/******************************************************************************
 * program:     rasimg library                                                *
 * function:    Miscellaneous operations with raster.                         *
 * modul:       imageop.cc                                                    *
 * licency:     GPL or LGPL                                                   *
 ******************************************************************************/
#include <stddef.h>
#include <stdio.h>

#include <stdlib.h>
#include <ctype.h>
#include <string.h>

#include "typedfs.h"
#include "image.h"
#include "ras_prot.h"
#include "imageop.h"
#include "img_tool.h"

#include "matrix.h"


//#include <alloc.h>

#ifndef MIN
 #define MIN(a,b) ((a)<(b) ? (a) : (b))
#endif


inline int32_t sqr(const int32_t val)
{
  return val*val;
}

inline float sqr(const float val)
{
  return val*val;
}



////////////////////////////////////////////////////////////////////

Image PeelPlane(const Image Img, int PlaneNum)
{
Raster2DAbstract *Raster=NULL;
Image NewImg;

 if(Img.Raster!=NULL)
  {
  Raster = CreateRaster2D(Img.Raster->Size1D,Img.Raster->Size2D,1);
  if(Raster)
    for(unsigned y=0; y<Img.Raster->Size2D; y++)	// main row loop
      {
      Img.Raster->GetRowRaster(y)->Peel1Bit(Raster->GetRow(y),PlaneNum);
      }
  }

  if(Raster!=NULL) NewImg.AttachRaster(Raster);
  if(Img.Palette!=NULL) NewImg.AttachPalette(Img.Palette);

return(NewImg);
}



Image Peel8bit(const Image Img, int PlaneNum8)
{
Raster2DAbstract *Raster=NULL;
Image NewImg;

 if(Img.Raster!=NULL)
  {
  PlaneNum8 *= 8;
  Raster = CreateRaster2D(Img.Raster->Size1D,Img.Raster->Size2D,8);

  if(Raster)
    for(unsigned y=0; y<Img.Raster->Size2D; y++)	// main row loop
      {
      Img.Raster->GetRowRaster(y)->Peel8Bit((unsigned char*)Raster->GetRow(y),PlaneNum8);
      }
  }

  if(Raster!=NULL) NewImg.AttachRaster(Raster);
  if(Img.Palette!=NULL) NewImg.AttachPalette(Img.Palette);

return(NewImg);
}


/** Peels one bit plane.
 * @param[in]	obr	Source image with multiple planes.
 * @param[out]	PlaneX	Preallocated 1 plane image.
 * @param[in]	PlaneNum	Number of a plane. */
void PeelPlane(Raster2DAbstract & obr, Raster2DAbstract & PlaneX, int PlaneNum)
{
 if(PlaneX.GetPlanes()!=1 || PlaneX.Size1D!=obr.Size1D) return;

 if(PlaneNum<0 || PlaneNum>=obr.GetPlanes() || obr.Data2D==NULL) return;
 //PlaneX.Create(obr.SizeX(),obr.SizeY(),1 | NoFill);

 unsigned MaxY = (obr.Size2D<PlaneX.Size2D)?obr.Size2D:PlaneX.Size2D;
 while(MaxY-- > 0)	// main row loop
 {
   obr.GetRowRaster(MaxY)->Peel1Bit(PlaneX.GetRow(MaxY),PlaneNum);
 }
}



void JoinPlane(const Raster2DAbstract & PlaneX, Raster2DAbstract & Obr, uint8_t PlaneNum)
{
 if(((PlaneNum<0)||(PlaneNum>=Obr.GetPlanes()))||(PlaneX.Data2D==NULL)) return;
 if(Obr.Data2D==NULL)
 {/*	//if no raster available, create it and hope that PlaneNum is a maximal plane
    *Obr=CreateRaster2D(PlaneX.Size1D,PlaneX.Size2D,PlaneNum);
    if(Obr.Data2D==NULL)*/
   return;
 }

 unsigned MaxY = (Obr.Size2D<PlaneX.Size2D) ? Obr.Size2D : PlaneX.Size2D;
 while(MaxY-- > 0)	// Main row loop
 {
    Obr.GetRowRaster(MaxY)->Join1Bit(PlaneX.GetRow(MaxY), PlaneNum);
 }
}


////////////////////////////////////////////////////////////////////


/// Filter image using a mask with given size.
Image Mask(Image *pIm_in, int size)
{
Image im_out, *pImOut=&im_out;
uint16_t k2;
uint32_t k;
uint16_t i,j;
uint16_t x,y;
int y1,y2,x1,x2;
uint32_t aLl;
uint32_t *ptrLine;
uint8_t *Sptr;
uint32_t *SptrAl,*SptrAr;
Raster1D_8Bit ptrLineB;
uint8_t color,clr;

 if(size==0)  return im_out;

 while(pIm_in!=NULL)
 {
   color = 1;
   if(pIm_in->ImageType()==ImageTrueColor)
   {
     color = 3;
     //  Obr2.typ:='C';
     //if AlineProc<>nil then AlineProc^.pases := AlineProc^.pases*3;
   }

   if(pImOut->Raster != NULL)
   {
     pImOut->Next = new Image;
     pImOut = pImOut->Next;
   }
   pImOut->Create(pIm_in->Raster->GetSize1D(), pIm_in->Raster->Size2D, pIm_in->Raster->GetPlanes());
   if(pImOut->Raster==NULL || pImOut->Raster->Data1D==NULL) break;

   //if( AlineProc!=NULL) AlineProc->InitPassing(Obr1.y,'Filtering by MASK');

   ptrLine = (uint32_t *)malloc(4*pIm_in->Raster->GetSize1D());
   ptrLineB.Allocate1D(pIm_in->Raster->GetSize1D());

   k = size*size;
   k2 = k/2;
   x = size/2 + 1;
   y = size/2 + 1;

   for(clr=0;clr<color;clr++)
   {
     y1 = x + 1;
     if(size<=2) y1--;		//pocet pridanych radek

     memset(ptrLine,0,4*pIm_in->Raster->GetSize1D());
     if(color==1)
        ptrLineB.Set(*pIm_in->Raster->GetRowRaster(0));
     else
       pIm_in->Raster->GetRowRaster(0)->Peel8Bit((uint8_t *)ptrLineB.Data1D,8*clr);


     for(j=0; j<y1; j++) 
       AddLu32u8(ptrLineB.GetSize1D(),ptrLine,(uint8_t *)ptrLineB.Data1D); /* ukradek=y1*ukradekW*/

     for(y1=1;y1<=(int)(size+1)/2 - 2;y1++)
     {		// pridani pripravnych radku}
       if(color==1)
         ptrLineB.Set(*pIm_in->Raster->GetRowRaster(y1));
       else
         pIm_in->Raster->GetRowRaster(y1)->Peel8Bit((uint8_t *)ptrLineB.Data1D,8*clr);
       AddLu32u8(ptrLineB.GetSize1D(),ptrLine,(uint8_t *)ptrLineB.Data1D);
     }

     for(i=0;i<pImOut->Raster->Size2D;i++)	// main filtration loop
     {
       y1 = (int)(i)-y;
       y2 = y1+size;
       if(y1<0) y1=0;
       if(y2>=pImOut->Raster->Size2D) y2 = pImOut->Raster->Size2D - 1;

      if(color==1)
        ptrLineB.Set(*pIm_in->Raster->GetRowRaster(y1));
     else
        pIm_in->Raster->GetRowRaster(y1)->Peel8Bit((uint8_t *)ptrLineB.Data1D,8*clr);
      SubLu32u8(ptrLineB.GetSize1D(),ptrLine,(uint8_t *)ptrLineB.Data1D); //vyjmuti posledniho radku y1

     if(color==1)
       ptrLineB.Set(*pIm_in->Raster->GetRowRaster(y2));
     else
       pIm_in->Raster->GetRowRaster(y2)->Peel8Bit((uint8_t *)ptrLineB.Data1D,8*clr);
     AddLu32u8(ptrLineB.GetSize1D(),ptrLine,(uint8_t *)ptrLineB.Data1D); //pridani noveho radku y2
		//zpracovani obsahu radku
     x1 = y + 1;
     if(size<=2) x1--;		// pocet pridanych sloupcu
     aLl = x1 * ptrLine[0];
     for(x2=1;x2<=(int)(size+1)/2 - 2;x2++) aLl += ptrLine[x2];
 
     Sptr = (uint8_t *)ptrLineB.Data1D;		//vypocet masky v radku
     x1 = -(int)(x);
     x2 = MIN(x1+size,ptrLineB.GetSize1D()-1);
     SptrAl = ptrLine;
     SptrAr = ptrLine + x2;
     for(j=0;j<ptrLineB.GetSize1D();j++)
     {
       aLl += *SptrAr;
       aLl -= *SptrAl;

       *Sptr = (aLl + k2) / k;

       if(x1>=0)	SptrAl++;
       if(x2<ptrLineB.GetSize1D()-1) SptrAr++;
       Sptr++;
       x1++;
       x2++;
     }

     if(color==1)		// data storage
       ptrLineB.Get(*pImOut->Raster->GetRowRaster(i));
     else
       pImOut->Raster->GetRowRaster(i)->Join8Bit((uint8_t *)ptrLineB.Data1D,8*clr);

     //if AlineProc<>nil then AlineProc^.NextLine;	{Dialog s uzivatelem}
    }
   }
   ptrLineB.Erase();
   free(ptrLine);  ptrLine=NULL;

   pIm_in = pIm_in->Next;
 }

return(im_out);
}


///////////////////////////////////////////////////////////////////////////

extern "C" {
typedef struct
{
  uint16_t MaxSz;
  uint16_t Size;
  uint16_t data[];
} TMed;

#ifdef _DEBUG
void PrintElement(const TMed *pMed)
{
  putchar('[');
  for(unsigned i=0; i<pMed->Size; i++)
  {
    if(i>0) putchar(',');
    printf("%u",pMed->data[i]);
  }
  putchar(']');
}
#endif


void InsertElement(uint16_t element, TMed *pMed)
{
unsigned i;
uint16_t *pData;
 if(pMed!=NULL && pMed->Size<pMed->MaxSz)
 {
   pData = pMed->data;
   i = ++pMed->Size;

   while(i>0)
   {
     i--;
     if(i==0)
     {
       pData[0] = element;
       break;
     }
     if(pData[i-1]<=element)
     {
       pData[i] = element;
       break;
     }
     pData[i] = pData[i-1];
   }
 }
}


void RemoveElement(uint16_t element, TMed *pMed)
{
uint16_t i;
uint16_t *pData;
 if(pMed)
 {
   pData = pMed->data;
   for(i=0; i<pMed->Size; i++)
   {
      if(pData[i]==element)
      {
        i++;
        while(i<pMed->Size)
        {
          pData[i-1] = pData[i];
          i++;
        }
        pMed->Size--;
        break;
      }
    }
   }
 }
}


void RemoveBlock(const TMed *pRmMed, TMed *pMed)
{
uint16_t i,j, r;
uint16_t *pData;
const uint16_t *pRmData;
 if(pMed==NULL || pRmMed==NULL) return;

 pData = &(pMed->data[0]);
 pRmData = &(pRmMed->data[0]);

 j = 0;
 r = 0;
 for(i=0; i<pMed->Size; i++)
 {
   if(r<pRmMed->Size && pData[i]==pRmData[r])
   {
     r++;
   }
   else
   {
     pData[j] = pData[i];
     j++;
   }
 }
 pMed->Size = j;
}


void InsertBlock(const TMed *pAddMed, TMed *pMed)
{
uint16_t i,j, a;
uint16_t *pData;
const uint16_t *pAddData;
 if(pMed==NULL || pAddMed==NULL) return;

 pData = &(pMed->data[0]);
 pAddData = &(pAddMed->data[0]);

 j = pMed->Size;
 a = pAddMed->Size;
 i = j + a;
 if(i>pMed->MaxSz) i=pMed->MaxSz;
 pMed->Size = i;

 while(i>0)
 {
   i--;
   if(j==0)
   {
     pData[i] = pAddData[--a];
     continue;
   }
   if(a==0)
   {
     //pData[i] = pData[--j];
     //continue;
     break;	// all data needad ara already here; i=j-1
   }

   if(pData[j-1] >= pAddData[a-1])
     pData[i] = pData[--j];
   else
     pData[i] = pAddData[--a];
 }
}


inline unsigned FitInterval(int x, int MaxX)
{
  if(x<=0) return 0;
  if(x>=MaxX) return MaxX-1;
  return (unsigned)x;
}


Image Median(const Image &im_in, int n)
{
int i,j;
int x,x1,y,y1;
uint16_t q;
uint8_t Color;
TMed *MediR;
TMed *MediG=NULL, *MediB=NULL;
TMed *Lmed;
TMed *LmedG=NULL, *LmedB=NULL;
Image im_out;
RGBQuad RGB;
const Raster1DAbstract *pRowLine;

 if(n==0 || im_in.Raster==NULL) return im_out;

 Color = 1;
 if(im_in.ImageType()==ImageTrueColor)
 {
   Color = 3;
   //  Obr2.typ:='C';
   im_out.Create(im_in.Raster->GetSize1D(), im_in.Raster->Size2D, ImageTrueColor|(im_in.Raster->GetPlanes()/im_in.Raster->Channels()));
 }
 else
   im_out.Create(im_in.Raster->GetSize1D(),im_in.Raster->Size2D,im_in.Raster->GetPlanes());

 if(im_out.Raster==NULL || im_out.Raster->Data1D==NULL) return im_out;

 // if AlineProc<>nil then AlineProc^.InitPassing(Obr1.x,'Filtering by MEDIAN');

 q = n*n;
 MediR = (TMed*)malloc(2*q+offsetof(TMed,data));
 Lmed = (TMed*)malloc(2*n+offsetof(TMed,data));
 if(MediR==NULL || Lmed==NULL) goto Fail;
 MediR->MaxSz = q;  
 Lmed->MaxSz = n;
 if(Color>1)
 {
   MediG = (TMed*)malloc(2*q+offsetof(TMed,data));
   MediB = (TMed*)malloc(2*q+offsetof(TMed,data));
   if(MediG==NULL || MediB==NULL) goto Fail;
   MediG->MaxSz = MediB->MaxSz = q;
   LmedG = (TMed*)malloc(2*n+offsetof(TMed,data));
   LmedB = (TMed*)malloc(2*n+offsetof(TMed,data));
   if(LmedG==NULL || LmedB==NULL) goto Fail;
   LmedG->MaxSz = LmedB->MaxSz = n;
 }

 q /= 2;
 x1 = y1 =  n / 2;

 {
   for(i=0; i<im_in.Raster->Size1D; i++)		// main filter loop
   {
     MediR->Size = 0;
     if(Color>1) MediG->Size=MediB->Size = 0;
     for(x=i-x1; x<i-x1+n; x++)
     {
       const int xx = FitInterval(x,im_in.Raster->Size1D);
       if(Color<=1)
       {
	 for(y=-y1-1; y<=-y1+n-2; y++)			// add one row before mask that will be removed
         {
	   //printf("IE:%d,%d,%u",x,y,im_in.Raster->GetValue2D(x,y));
	   InsertElement(im_in.Raster->GetValue2D(xx,(y<0)?0:y),MediR);
         }
       }
       else
       {         
 	 for(y=-y1-1; y<=-y1+n-2; y++)
         {
           im_in.Raster->Get(xx, (y<0)?0:y, &RGB);
	   InsertElement(RGB.R,MediR);
           InsertElement(RGB.G,MediG);
           InsertElement(RGB.B,MediB);
         }
        }
     }

     //if(n>5)		// Plain unblocked algoritmh has no benefit in C++, removed.
     for(j=0; j<im_in.Raster->Size2D; j++)
     {
       y = j-y1-1;
       pRowLine = im_in.Raster->GetRowRaster((y<0)?0:y);  
       if(pRowLine==NULL) continue;
       Lmed->Size = 0;

       if(Color<=1)
       {         
 	 for(x=i-x1; x<i-x1+n; x++) InsertElement(pRowLine->GetValue1D(FitInterval(x,im_in.Raster->Size1D)),Lmed);
       }
       else
       {
         LmedB->Size = LmedG->Size = 0;
	 for(x=i-x1; x<i-x1+n; x++) 
         {
           pRowLine->Get(FitInterval(x,im_in.Raster->Size1D), &RGB);
           InsertElement(RGB.R, Lmed);
           InsertElement(RGB.G, LmedG);
           InsertElement(RGB.B, LmedB);
         }
         RemoveBlock(LmedG,MediG);
         RemoveBlock(LmedB,MediB);
       }
        
       RemoveBlock(Lmed,MediR);

       y += n;
       pRowLine = im_in.Raster->GetRowRaster((y>=im_in.Raster->Size2D)?(im_in.Raster->Size2D-1):y);
       if(pRowLine==NULL) continue;
       Lmed->Size = 0;
       if(Color<=1)
       {
         for(x=i-x1; x<i-x1+n; x++)
            InsertElement(pRowLine->GetValue1D(FitInterval(x,im_in.Raster->Size1D)),Lmed);
       }
       else
       {
         LmedB->Size = LmedG->Size = 0;
         for(x=i-x1; x<i-x1+n; x++)
         {
           pRowLine->Get(FitInterval(x,im_in.Raster->Size1D), &RGB);
           InsertElement(RGB.R,Lmed);
           InsertElement(RGB.G,LmedG);
           InsertElement(RGB.B,LmedB);
         }
         InsertBlock(LmedG,MediG);
         InsertBlock(LmedB,MediB);
       }

        //PrintElement(Lmed);PrintElement(MediR);
       InsertBlock(Lmed,MediR);
        //PrintElement(MediR);

       if(Color<=1)
	  im_out.Raster->SetValue2D(i,j,MediR->data[q]);		// data store
       else
       {
         RGB.R = MediR->data[q];
         RGB.G = MediG->data[q];
         RGB.B = MediB->data[q];
         im_out.Raster->Set(i,j,&RGB);
       }
      }
      //if AlineProc<>nil then AlineProc^.NextLine;
   }
 }

Fail:
 if(Lmed) free(Lmed);
 if(LmedG) free(LmedG);
 if(LmedB) free(LmedB);
 if(MediG) free(MediR);
 if(MediG) free(MediG);
 if(MediB) free(MediB);

 return im_out;
}


///////////////////////////////////////////////////////////////////////////

Image Threshold(const Image *pImg, unsigned ThrValue)
{
Image NewImg, *pNewImg=&NewImg;
unsigned x, y;

 while(pImg != NULL)
 {
   if(pImg->Raster!=NULL && pImg->Raster->Size1D>0 && pImg->Raster->Size2D!=0)
   {
     if(pNewImg->Raster!=NULL)
     {
       pNewImg->Next = new Image;
       pNewImg = pNewImg->Next;
     }
     pNewImg->AttachRaster(CreateRaster2D(pImg->Raster->Size1D,pImg->Raster->Size2D,1));     
     if(pNewImg->Raster)
       for(y=0; y<pImg->Raster->Size2D; y++)		// main row loop
       {
         Raster1DAbstract * const Rsrc = pImg->Raster->GetRowRaster(y);
         Raster1DAbstract * const Rthr = pNewImg->Raster->GetRowRaster(y);

         for(x=0; x<pImg->Raster->Size1D; x++)
         {
           if(Rsrc->GetValue1D(x)>ThrValue)
              Rthr->SetValue1D(x,1);
           else
              Rthr->SetValue1D(x,0);
         }
     }
   }   
   pImg = pImg->Next;
 }
return(NewImg);
}


///////////////////////////////////////////////////////


void Convert2Gray(Image *pI)
{
Raster2DAbstract *Raster = NULL;

  while(pI!=NULL)
  {
    switch(pI->ImageType())
    {
      case ImageNone:
      case ImageGray:
      case ImagePaletteF:
      case ImageReal:
      case ImageSigned:		break;
      
      case ImagePalette:	 //2-palette,
		Raster = CreateRaster2D(pI->Raster->Size1D,pI->Raster->Size2D,8);
                if(Raster)
                {
                  for(unsigned y=0; y<Raster->Size2D; y++)		// main row loop
                  {
                    Raster1DAbstract *RSrc = pI->Raster->GetRowRaster(y);
                    Raster1DAbstract *RDest = Raster->GetRowRaster(y);
                    for(unsigned x=0; x<Raster->Size1D; x++)
                    {
                      uint32_t val = pI->Palette->GetValue1D(RSrc->GetValue1D(x));
                      val = ((val&0xFF) + ((val>>8)&0xFF) + ((val>>16)&0xFF)) / 3;
	              RDest->SetValue1D(x,val);
                    }
                  }
                  pI->AttachRaster(Raster);
                  pI->AttachPalette(NULL);
                  Raster = NULL;
                }
                break;

      case ImageTrueColor:
                Raster = CreateRaster2D(pI->Raster->Size1D,pI->Raster->Size2D,8);
                if(Raster)
                {
                  for(unsigned y=0; y<Raster->Size2D; y++)		// main row loop
                  {
                    Raster1DAbstract *RSrc = pI->Raster->GetRowRaster(y);
                    Raster1DAbstract *RDest = Raster->GetRowRaster(y);
                    for(unsigned x=0; x<Raster->Size1D; x++)
                    {
                      uint32_t val = RSrc->GetValue1D(x);
                      val = ((val&0xFF) + ((val>>8)&0xFF) + ((val>>16)&0xFF)) / 3;
	              RDest->SetValue1D(x,val);
                    }
                  }
                  pI->AttachRaster(Raster); Raster=NULL;  // Ussage count is 1; piinter could be discarded
                  pI->AttachPalette(NULL);                  
                }
                break;
		
    }
    pI = pI->Next;
  }
}


/// Canny edge detection for one pixel.
/// param[in]	BoxG	Surrounding matrix 3x3 with sqrt(sqr(ux)+sqr(uy)).
/// param[in]	ux32	x difference in a given pixel, packet as offseted 0 ~ 0x80000000.
/// param[in]	uy32	y difference in a given pixel, packet as offseted 0 ~ 0x80000000.
static inline bool Edge(const uint32_t * const *BoxG, const uint32_t ux32, const uint32_t uy32)
{
 if(ux32==0x80000000 && uy32==0x80000000) return false;

 const float fux = (int32_t)(ux32-0x80000000);
 const float fuy = (int32_t)(uy32-0x80000000);
 const float gx = fabs(fux * BoxG[1][1]);
 const float gy = fabs(fuy * BoxG[1][1]);

 if((ux32>=0x80000000) ^ (uy32<0x80000000))		// sign(GausDiffX)==sign(GausDiffX)
 {
   const float uz = fux - fuy;
   if(fabs(fux) < fabs(fuy))	//hrana ve smeru X
   {
     if((gy < fabs(fux*BoxG[2][2] - uz*BoxG[1][2])) ||
	(gy <= fabs(fux*BoxG[0][0] - uz*BoxG[1][0]))) return true;
   }
   else
   {				//hrana ve smeru Y
     if((gx < fabs(fuy*BoxG[2][2] + uz*BoxG[2][1])) ||
	(gx <= fabs(fuy*BoxG[0][0] + uz*BoxG[0][1]))) return true;
   }
 }
 else							//  {sign(GausDiffX)!=sign(GausDiffX)}
 {
   const float ut = fux + fuy;
   if(fabs(fux) < fabs(fuy))	// hrana ve smeru X
   {
     if((gy < fabs(fux*BoxG[2][0] - ut*BoxG[1][0])) ||
	(gy <= fabs(fux*BoxG[0][2] - ut*BoxG[1][2]))) return true;
   }
   else				// {hrana ve smeru Y}
   {
     if((gx < fabs(fuy*BoxG[2][0] - ut*BoxG[2][1])) ||
	(gx <= fabs(fuy*BoxG[0][2] - ut*BoxG[0][1]))) return true;
   }
 }
 return false;
}


/////////////////////////////////////////////////////////////////////////////////////////////


// This shifts lines line[0]=line[1]; line[1]=line[2]; ..... line[n-2]=line[n-1]; line[n-1]=line[0];
void RollLine(Raster2DAbstract * const Ras)
{
  if(Ras==NULL || Ras->Data2D==NULL || Ras->Size2D<=1) return;
  void * const Ptr0 = Ras->Data2D[0];
  for(unsigned y=1; y<Ras->Size2D; y++)
  {
    Ras->Data2D[y-1] = Ras->Data2D[y];
  }
  Ras->Data2D[Ras->Size2D-1] = Ptr0;
}


#define GausDiffY \
   {I = SizeY-GausN-1;		\
   void * const PLine = LB2.Data2D[I];	\
   LB2.Data2D[I] = LB1.Data2D[I-1];	\
   LB1.Data2D[I-1] = PLine;	\
   DiffVer1_u32((uint32_t*)PLine, (const uint32_t*)LB2.Data2D[I-1], (const uint32_t*)LB1.Data2D[I], SizeX); }


void GausFunc(Raster2DAbstract *obr, unsigned char FuncC, unsigned char GausN)
{
uint8_t DiffN;
int  I,YI,YO,SizeY;
Raster2D_32Bit LB1, LB2;
Raster2D_32Bit R3L;
//float CannyScale = 1;
Raster1D_8Bit PL8;
uint8_t Clr;

  DiffN=0;
  switch(FuncC)
  {
    case gfGaus:    break;
    case gfGausDiffX:
    case gfGausDiffY:
    case gfGausDiff: 
    case gfCanny:
    case gfCannyFast:
    case gfGausLap  :DiffN=1;
		     break;
  }
  //if(Mult==0) Mult=0x10000;
  const int SizeX = obr->Size1D;
  if(SizeX<10 || obr->Size2D<10) return;

  if((FuncC!=gfCanny)&&(FuncC!=gfCannyFast))
    SizeY=1+GausN+DiffN;
  else
  {
    SizeY = 2+GausN+DiffN;
    R3L.Allocate2D(SizeX,3);
    for(I=0; I<R3L.Size2D; I++)
    {
      memset(R3L.GetRow(I),0,SizeX*sizeof(float));
    }
    //CannyScale=Mult;
    //CannyScale=4/CannyScale;
  }

  LB1.Allocate2D(SizeX,SizeY);
  LB2.Allocate2D(SizeX,SizeY);

  const uint8_t Colors = obr->Channels();
  if(Colors>0) PL8.Allocate1D(SizeX);

  for(Clr=0; Clr<Colors; Clr++)
  {    
    YI = -(SizeY-1);
    YO = -(2*SizeY-2);
    while(YO < (int)obr->Size2D)
    {    
      RollLine(&LB1);
      RollLine(&LB2);

      Raster1DAbstract *PL = LB1.GetRowRaster(SizeY-1);
	// Duplicate first/last row only - no reflection.
    
      if(YI<0) I=0;
      else if(YI>=obr->Size2D) I=obr->Size2D-1;
      else I=YI;      

      if(Colors==1)
         PL->Set(*obr->GetRowRaster(I));
      else
      {
        obr->GetRowRaster(I)->Peel8Bit((uint8_t *)PL8.Data1D,8*Clr);
        PL->Set(PL8);        
      }


/*	// Reflection - commented out
    if(YI<(int)obr->Size2D)
      PL->Set(*obr->GetRowRaster(labs(YI)));
    else
    {
      int YY = 2*obr->Size2D-YI-2;
      if(YY<0) YY=0;      
      PL->Set(*obr->GetRowRaster(YY));
    } */

    Gaus1D_u32((uint32_t*)PL->Data1D,SizeX,GausN);
    YI++;

    for(I=SizeY-1; I>=SizeY-GausN; I--)
    {
      void *PLine = LB2.Data2D[I];
      LB2.Data2D[I] = LB1.Data2D[I-1];
      LB1.Data2D[I-1] = PLine;
      GausVer1_u32((uint32_t*)PLine, (uint32_t*)LB1.Data2D[I], (uint32_t*)LB2.Data2D[I], (uint32_t*)LB2.Data2D[I-1], SizeX);
    }

    switch(FuncC)
    {
      case gfGaus:if(YO>=0)
                  {
                    PL = LB1.GetRowRaster(0);
                    if(Colors==1)
                      obr->GetRowRaster(YO)->Set(*PL);
                    else
                    {
                      PL8.Set(*PL);
                      obr->GetRowRaster(YO)->Join8Bit(PL8.Data1D,8*Clr);
                    }
                  }
                  break;

      case gfGausDiffX:if(YO>=0)
                  {
	            PL = LB1.GetRowRaster(0);
	            DiffHor1_u32((uint32_t*)PL->Data1D, (const uint32_t*)PL->Data1D, SizeX);
                    if(Colors==1)
	              obr->GetRowRaster(YO)->Set(*PL);
                    else
                    {
                      PL8.Set(*PL);
                      obr->GetRowRaster(YO)->Join8Bit(PL8.Data1D,8*Clr);
                    }
                  }
                  break;

      case gfGausDiffY:
	          GausDiffY
                  if(YO>=0)
                  {
	            PL = LB1.GetRowRaster(0);
                    if(Colors==1)
	                obr->GetRowRaster(YO)->Set(*PL);
                    else
                    {
                      PL8.Set(*PL);
                      obr->GetRowRaster(YO)->Join8Bit(PL8.Data1D,8*Clr);
                    }
                  }
                  break;


      case gfGausDiff:
	         GausDiffY
	         if(YO>=0)
	         {
	           I = SizeY - GausN - 2;	           
	           DiffHor1_u32((uint32_t*)LB2.GetRow(I), (uint32_t*)LB2.GetRow(I+1), SizeX);
	           AbsDiff_u32((uint32_t*)LB1.GetRow(I), (uint32_t*)LB2.GetRow(I), SizeX);
	           PL = LB1.GetRowRaster(0);
                   if(Colors==1)
	               obr->GetRowRaster(YO)->Set(*PL);
                   else
                   {
                     PL8.Set(*PL);
                     obr->GetRowRaster(YO)->Join8Bit(PL8.Data1D,8*Clr);
                   }
	         }
                 break;

      case gfCanny:          
	      GausDiffY;
	      I = SizeY-GausN-2;
              DiffHor1_u32((uint32_t*)LB2.GetRow(I), (const uint32_t*)LB2.GetRow(I+1), SizeX);
              {
                uint32_t MaxDiff = 0;
                {
	          uint32_t * const PL1_32 = (uint32_t*)LB1.GetRow(I);
	          uint32_t * const PL2_32 = (uint32_t*)LB2.GetRow(I);
                  uint32_t * const PRL = (uint32_t*)R3L.GetRow(0);
	          for(I=0; I<SizeX; I++)
	          {
	            PRL[I] = sqrt((float)(sqr(((float)((int)(PL1_32[I]-0x80000000)))) + sqr((float)((int)(PL2_32[I]-0x80000000)))));
                    if(PRL[I]>MaxDiff) MaxDiff=PRL[I];
	          }
                }

                RollLine(&R3L);	        

	        if(YO>=0)
	        {
                  uint8_t * const pPL8 = (uint8_t *)PL8.Data1D;
	          uint32_t * const PL1_32 = (uint32_t*)LB1.GetRow(0);
	          uint32_t * const PL2_32 = (uint32_t*)LB2.GetRow(0);
                  uint32_t * R3L_32[3] = {(uint32_t*)R3L.GetRow(0), (uint32_t*)R3L.GetRow(1), (uint32_t*)R3L.GetRow(2)};
                  MaxDiff = 0xFFFFFFFF / MaxDiff;
                  if(MaxDiff==0) MaxDiff=1;
	          for(I=1; I<=SizeX-2; I++)
	          {
	            if(!Edge(R3L_32, PL1_32[I], PL2_32[I]))
	              pPL8[I] = 0;
	            else
	            {
	              uint32_t R1 = (MaxDiff*R3L_32[1][1]) >> 24;
	              pPL8[I] = (R1>255) ? 255 : R1;
	            }
                    R3L_32[0]++; R3L_32[1]++; R3L_32[2]++;
	          }
	          pPL8[0] = pPL8[SizeX-1] = 0;                

                  if(Colors==1)
	            obr->GetRowRaster(YO)->Set(PL8);
                  else
                    obr->GetRowRaster(YO)->Join8Bit(PL8.Data1D,8*Clr);
                }
	      }
              break;

      case gfCannyFast:
	GausDiffY
	I = SizeY-GausN-2;        
	DiffHor1_u32((uint32_t*)LB2.GetRow(I), (const uint32_t*)LB2.GetRow(I+1), SizeX);	
	AbsDiffCopy_u32((uint32_t*)R3L.GetRow(0), (uint32_t*)LB1.GetRow(1), (uint32_t*)LB2.GetRow(1), SizeX);
        RollLine(&R3L);
	if(YO>=0)
	{
//	  PL1 = LB1.Data2D[0];
//	  PL2 = LB2.Data2D[0];

//	  PL =PL1;		{ Funkce EdgeLine v Assembleru }
	  //EdgeLine(LB1.Data2D[0],pointer(@R3L),PL1,PL2,SizeX);


          PL = LB1.GetRowRaster(0);
          if(Colors==1)
	      obr->GetRowRaster(YO)->Set(*PL);
          else
          {
            PL8.Set(*PL);
            obr->GetRowRaster(YO)->Join8Bit(PL8.Data1D,8*Clr);
          }
	}
        break;

/*
      case gfGausLap:{
	I=SizeY-GausN-1;
	PL =LB2.Data^[I];
	PL2=LB1.Data^[I-1];
	LB1.Data^[I-1]=PL;
	LB2.Data^[I]=PL2;
	DiffVer2(PL,LB2.Data^[I-1],PL2,LB1.Data^[I],SizeX);
	if YO>=0 
	{
	  I=SizeY-GausN-2;
	  PL =LB2.Data^[I];
	  PL1=LB2.Data^[I+1];
	  DiffHor2(PL,PL1,SizeX);
	  PL =LB1.Data^[I];
	  PL1=LB2.Data^[I];
	  AddHor(PL,PL1,SizeX);
	  WriteLine(obr,YO,PL,Mult,Offs);
	}
      }
*/
      }
      //if AlineProc!=nil  AlineProc^.NextLine;
      YO++;
    }
  }

  R3L.Erase();
 
  LB1.Erase();
  LB2.Erase();
}


#if 0
static const uint8_t removal[] = {
			     0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80,	// jeden soused
      			     0x03,0x06,0x0C,0x18,0x30,0x60,0xC0,0x81,	// dvojice sousedu
			     0x05,0x14,0x50,0x41,			// dvojice s dirou
      			     0x07,0x0E,0x1C,0x38,0x70,0xE0,0xC1,0x83,	// trojice
                             0x0D,0x34,0xD0,0x43,			// trojice s dirou
                             0x85,0x16,0x58,0x61,
                             0x15,0x54,0x51,0x45,			// trojice se 2 dirama
      			     0x0F,0x1E,0x3C,0x78,0xF0,0xE1,0xC3,0x87,	// ctverice
                             0x1D,0x74,0xD1,0x47,			// ctverice s dirou
			     0x8D,0x36,0xD8,0x63,
			     0xC5,0x17,0x5C,0x71,
                             0x35,0xD4,0x53,0x4D,			// ctverice se 2 dirama
                             0x95,0x56,0x59,0x65,
                             0x1F,0x3E,0x7C,0xF8,0xF1,0xE3,0xC7,0x8F,	// petice
                             0x3D,0xF4,0xD3,0x4F,			// petice s dirou
                             0x9D,0x76,0xD9,0x67,
			     0xCD,0x37,0xDC,0x73,
			     0xE5,0x97,0x5E,0x79,
			     0x3F,0x7E,0xFC,0xF9,0xF3,0xE7,0xCF,0x9F,	// sestice
			     0x00};					// nema okoli

void CannyCleanup(Raster2DAbstract *obr)
{
Raster2D_8Bit R2Cache;
Raster1D_8Bit Rout;
unsigned i;

  if(obr==NULL || obr->Size1D<=1 || obr->Size2D<=1) return;
  const unsigned SizeX = obr->Size1D;
  R2Cache.Allocate2D(SizeX,3);
  Rout.Allocate1D(SizeX);

  const uint8_t Colors = obr->Channels();
  
  for(uint8_t Clr=0; Clr<Colors; Clr++)
  {
    if(Colors==1)
        R2Cache.GetRowRaster(1)->Set(*obr->GetRowRaster(0));
    else
        obr->GetRowRaster(0)->Peel8Bit((uint8_t *)R2Cache.GetRow(1),8*Clr);
    memcpy(R2Cache.GetRow(2), R2Cache.GetRow(1), SizeX);

    for(unsigned y=0; y<obr->Size2D; y++)
    {
      unsigned Yp1 = y + 1;
      if(Yp1 >= obr->Size2D) Yp1=obr->Size2D-1;
      RollLine(&R2Cache);

      if(Colors==1)
          R2Cache.GetRowRaster(2)->Set(*obr->GetRowRaster(Yp1));
      else
          obr->GetRowRaster(Yp1)->Peel8Bit((uint8_t *)R2Cache.GetRow(2),8*Clr);
      
      Rout.SetValue1D(0, R2Cache.GetValue2D(0,1));
      for(unsigned x=0; x<SizeX; x++)
      {
        uint8_t a = R2Cache.GetValue2D(x,1);
	uint8_t m2=0;
        uint8_t Mask2=0;

        if(R2Cache.GetValue2D(x+1,1)!=0) {Mask2|= 1;m2++;}
        if(R2Cache.GetValue2D(x+1,0)!=0) {Mask2|= 2;m2++;}
        if(R2Cache.GetValue2D(x  ,0)!=0) {Mask2|= 4;m2++;}
        if(R2Cache.GetValue2D(x-1,0)!=0) {Mask2|= 8;m2++;}
        if(R2Cache.GetValue2D(x-1,1)!=0) {Mask2|=16;m2++;}
        if(R2Cache.GetValue2D(x-1,2)!=0) {Mask2|=32;m2++;}
        if(R2Cache.GetValue2D(x  ,2)!=0) {Mask2|=64;m2++;}
        if(R2Cache.GetValue2D(x+1,2)!=0) {Mask2|=128;m2++;}

        if(m2==1)    		//osamocene konce hran
        {
          if(GausN<15)
          {
            //if (GausN<6) continue; //nic uz neporoste
            if(a==0 && a2>0)	//pro male GausN potrebuji >=2 sousedy}
            {
              if((Mask2&0x55) != 0) continue;	//expanduj pouze do rohu

		                j=0;
                                if pixel(tempP,x+1,y  )!=0  j++;
		                if pixel(tempP,x+1,y-1)!=0  j++;
             	                if pixel(tempP,x  ,y-1)!=0  j++;
		                if pixel(tempP,x-1,y-1)!=0  j++;
             	                If pixel(tempP,x-1,y  )!=0  j++;
		                if pixel(tempP,x-1,y+1)!=0  j++;
             	                if pixel(tempP,x  ,y+1)!=0  j++;
             	                if pixel(tempP,x+1,y+1)!=0  j++;
                                if j<2  Continue; //jinak zustanou bez upravy
                                };
                        };

                     if a2>a 
                     	    {
			    SetPixel(im_out,x,y,a2);	//prodluz konec primky
                            Continue;
                            }
                     };		//jinak je vyradi nasledujici operace}


/*
             a=pixel(im_out,x,y);
             a2=pixel(tempP,x,y);

	     mask=$FF;mask2=$0;
             m2=0;
             if pixel(im_out,x+1,y  )!=0  {;mask2=Mask2 or  1;inc(m2);};
             if pixel(im_out,x+1,y-1)!=0  {;mask2=Mask2 or  2;inc(m2);};
             if pixel(im_out,x  ,y-1)!=0  {;mask2=Mask2 or  4;inc(m2);};
             if pixel(im_out,x-1,y-1)!=0  {;mask2=Mask2 or  8;inc(m2);};
             if pixel(im_out,x-1,y  )!=0  {;mask2=Mask2 or 16;inc(m2);};
             if pixel(im_out,x-1,y+1)!=0  {;mask2=Mask2 or 32;inc(m2);};
             if pixel(im_out,x  ,y+1)!=0  {;mask2=Mask2 or 64;inc(m2);};
             if pixel(im_out,x+1,y+1)!=0  {;mask2=Mask2 or 128;inc(m2);};

             if m2=1    		{osamocene konce hran}
             	     {
                     if (GausN<15) 
                     	{
                        if (GausN<6)  Continue; {nic uz neporoste}
                        if (a=0)and(a2>0) {pro male GausN potrebuji >=2 sousedy}
                        	{
                           	if Mask2 and $55 !=0  continue; {expanduj pouze do rohu}

		                j=0;
                                if pixel(tempP,x+1,y  )!=0  inc(j);
		                if pixel(tempP,x+1,y-1)!=0  inc(j);
             	                if pixel(tempP,x  ,y-1)!=0  inc(j);
		                if pixel(tempP,x-1,y-1)!=0  inc(j);
             	                If pixel(tempP,x-1,y  )!=0  inc(j);
		                if pixel(tempP,x-1,y+1)!=0  inc(j);
             	                if pixel(tempP,x  ,y+1)!=0  inc(j);
             	                if pixel(tempP,x+1,y+1)!=0  inc(j);
                                if j<2  Continue; {jinak zustanou bez upravy}
                                };
                        };

                     if a2>a 
                     	    {
			    SetPixel(im_out,x,y,a2);	{prodluz konec primky}
                            Continue;
                            }
                     };		{jinak je vyradi nasledujici operace}


             if a2=0 
             	{
                if (m2>=3)and(a=0) 		{korekce uzlu}
                   {			{male kolecko - vzdy vypln jeho stred}
                   if (Mask2 and $55 = $55) or	{Uzel S=1}
		       (((((Mask2 and $D6=$D6)and(pixel(im_out,x+2,y )!=0)) or {Uzel S=2}
                         (((Mask2 and $5B=$5B)and(pixel(im_out,x ,y+2)!=0)))) or
                        ((((Mask2 and $6D=$6D)and(pixel(im_out,x-2,y )!=0)) or
                         (((Mask2 and $B5=$B5)and(pixel(im_out,x ,y-2)!=0)))))))
                        	
                                {
                                if a<a2  a=a2;
                                if a=0  a=Pixel(im_out,x+1,y);
                                if a=0  a=1;
                                SetPixel(im_out,x,y,a);
                                Continue;
                                };
                   };

                if (mask2 in vyrazeni)	{vyjimani bodu}
		                  
                                  {
				  SetPixel(im_out,x,y,0);
                                  continue;
                                  };

                if a>1 		{drzeny bod- sniz duveru}
                	{
                        dec(a);
                        SetPixel(im_out,x,y,a);
                        };

                };


             if (a2!=0)and(m2>1) 
                {			{narustani hran, nebo}
                i=a2;			{totozna hrana = zvys pouze duveru}
                if a>i  i=a;
                SetPixel(im_out,x,y,i);

{		PutPixel(x,y,Yellow);{}
                Continue;
                };
*/

        Rout.SetValue1D(x, R2Cache.GetValue2D(x,1));
      }
      Rout.SetValue1D(SizeX-1, R2Cache.GetValue2D(SizeX-1,1));

      if(Colors==1)
          R2Cache.GetRowRaster(1)->Get(*obr->GetRowRaster(Yp1));
      else
          obr->GetRowRaster(Yp1)->Join8Bit((uint8_t *)R2Cache.GetRow(1),8*Clr);
    }
  }
}

#endif

